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Abstract 

Let f(x) be one of the usual elementary functions (exp, log, 
artan, sin, cosh, etc.), and let M(n) be the n-umber of single- 
precision operations required to multiply n-bit integers. We show 
that f(x) can be evaluated, -with relative error 0(2 ) , in 
0(M(n) log(n)) operations as n -» » , for any floating-point number 
X (-with an n-bit fraction) in a suitable finite interval. From the 
Schonhage-Strassen bound on M(n) , it follows that an n-bit approxi- 
mat ion to f(x) may be evaluated in 0(n log (n) log log (n)) 
operations. Special cases include the evaluation of constants such 
as jt , e , and e^ . The algorithms depend on the theory of elliptic 
integrals, using the arithmetic -geometric mean iteration and ascending 
Landen transformations. 
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1. Introduction . 

We consider the mjmber of operations required to evaluate the 
elementary functions exp(x) , log(x) -^ , artan(x) , sin(x) , etc., 
with relative error 0(2 ) , for x in some interval [a,b] , and 
large n . Here, [a^b] is a fixed, nontrivial interval on which 
the relevant elementary function is defined. The results hold for 
computations performed on a multitape Turing machine, but to simplify 
the exposition we assume that a standard serial computer with a 
random-access memory is used. 

Let M(n) be the number of operations required to multiply two 
integers in the range [0, 2 ) . We assume the number representation 
is such that addition can be performed in 0(M(n)) operations, and 
that M(n) satisfies the weak regularity condition 

M(an) < pM(n) , (1.1) 

for some a and p in (0,1) , and all sufficiently large n . 
Similar, but stronger, conditions are usually assiuned either explicitly 
[6] or implicitly [9]. Our assumptions are certainly valid if the 
Schonhage-Strassen method [9,11] is used to multiply n-bit integers 
(in the usual binary representation) in 0(n log(n) log log(n)) 
operations. 

The elementary function evaluations may be performed entirely in 
fixed point, using integer arithmetic and some implicit scaling scheme. 
However, it is more convenient to assume that floating-point computation 
is used. For example, a sign and magnitude representation could be 
used, with a fixed-length binaiy exponent and an n-bit binary fraction. , 
Our results are independent of the particular floating-point number 
system are used, so long as the following conditions are satisfied. 

1. Real numbers which are not too large or small can be approximated 
by floating-point numbers, with a relative error 0(2 ) . 

2. Floating-point addition and multiplication can be performed in 
0(M(n)) operations, with a relative error 0(2 ) in the result. 



W 



log(x) denotes the natural logarithm. 



5. The precision n is variable, and a fJ.oating -point naiaiber Mth 
precision n may be approximated, with relative error 0(2""^) 
and in 0(M(n)) operations, by a floating-point number with 
precision m , for any positive m < n . 

Throughout this paper, a floating-point number means a number in 
some representation satisfying conditions 1 to 5 above, not a single- 
precision number. We say that an operation is perfoimed with 
precision n if the result is obtained with a relative error 0(2 ) . 
It is assumed that the operands and result are approximated by floating- 
point numbers. 

The main result of this paper, established in Sections 6 and "J, 
is that all the usual elementary functions may be evaluated, with 
precision n , in 0(M(n) log(n)) operations. Note that 0(M(n)n) 
operations are required if the Taylor series for log(l+x) is summed 
in the obvious way. Our result improves the bound 0(M(n) log (n)) 
given in [3]> although the algorithms described there may be faster 
for small n . 

Preliminary results are given in Section 2 to 5. In Section 2 
we give, for completeness, the known result that division and 
extraction of square roots to precision n require 0(M(n)) operant ions. 
Section 5 deals briefl^r with methods for approximating simple zeros 
of nonlinear equations to precision n , and some results from the 
theory of elliptic integrals are summarized in Section h. Since our 
algorithm? for elementary functions require a knowledge of jt to 
precision n , we show, in Section ^, how this may be obtained in 
0(M(n) log(n)) operations. An amusing consequence of the results of 
Section 6 is that e may also be evaluated, to precision n , in 
0(M(n) log(n)) operations. 

From Theorem ^.1 of [3]y at least 0(M(n)) operations are 
required to evaluate exp(x) or sin(x) to precision n . It is 
plausible to conJect\ire that 0(M(n) log(n)) operations are necessary. 

Most of this paper is concerned with order of magnitude results, 
and multiplicative constajits are ignored. In Section 8, though, we 
give upper boiaids on the constants. From these bounds it is possible 
to estimate how large n needs to be before our algorithms are faster 
than the conventional ones. 
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After this report was written^ Bill Gosper drew my attention to 
the paper of Salamin [10], where an algorithm very similar to o\ir 
algorithm for Evaluating n is descrihed. A fast algorithm for 
evaluating log(x) was also found (independently) by Salarain 
(see [2]). 



2. Reciprocals and Square Roots . 

In this section we show that reciprocals and square roots of 
floating-point numbers may be evaluated, to precision n , in 0(M(n)) 
operations . To simplify the statement of the following lemma, we 
assume that M(x) = for all x < 1 . 

00 

Lemma 2.1 . If 7 e (0,1) , then E M(7'^n) = 0(M(n)) as n - » . 

Proof . If a and p are as in (1.1), there exists k such that 

00 , 00 ^ 

y^ <a , Thus, L M(7^n) <k S U{oPn) < kM(n)/(l-p) + 0(1) , 

j=0 j=0 

by repeated application of (1.1) . Since M(n) -» «> as n -» °° , the 

result follows. Q 

In the following lemma, we assume that l/c is in the allowable 
range for floating-point nimibers. Similar assumptions are implicit 
below. 

Lemma 2.2 . If c is a nonzero floating-point number, then l/c can 
be evaluated, to precision n , in 0(M(n)) operations. 

Proof . The Newton iteration 

x.._ = x.(2 -ex.) / (2.1) 

1+1 1^ 1' ' V / 

converges to l/c with order 2 . In fact, if x. = (l-e.)/c , 
substitution in (2.1) gives e. -, = e. . Thus, assimiing je^l < l/2 , 



we have |e.| < 2 for all i >0 , and x is a sufficiently 
good approximation to l/c if k > log^ n . This ass-umes that (2.1) 
is satisfied exactly, but it is easy to show that it is sufficient to 
use precision n at the last iteration (i = k-1) j precision 
slightly greater than n/2 for i = k-2 , etc. (Details, and more 
efficient methods, are given in [3jh].) Thus, the result follows from 
Lemma 2.1. Since x/y = x(l/y) , it is clear that floating-point 
division may also be done in 0(M(n)) operations. D 

1/2 
Lemma 2.3 * If '^ ^ ^ is a floating-point number, then c ' can be 

evaluated, to precision n , in 0(M(n)) operations. 

1/2 
Proof . If c = then c ' =0. Ifc^O, the proof is similar 

to that of Lemma 2.2, using the Newton iteration x. ^ = (x. + c/x.)/2 . Q 
Lemma 2.U . For £iny fixed k > , M(kn) = 0(M(n)) as n -♦ 0° . 

Proof . Since we can add integers less than 2 in 0(M(n)) 

operations, we can add integers less than 2 in 0(kM(n)) = 0(M(n)) 

kn 
operations . The multiplication of integers less than 2 can be 

2 n 

split into 0(k ) multiplications of integers less than 2 , and 

0(k ) additions, so can be done in 0(k'li(n)) = 0(M(n)) operations. Q 



5* Solution of Nonlinear Equations . 

In Section 6 we need to solve nonlinear equations to precision n . 
The following lemma is sufficient for this application. Stronger results 
are given in [5^^]» 

Lemma 3'1 » If the equation f (x) = c has a simple root ^ /= , 
f is Lipschitz continuous near ^ , and we can evaluate f (x) to 
precision n in 0(M(n)^(n)) operations, where 0(n) is a positive, 
monotonic increasing function, for x near Q , then Q can be 
evaluated to precision n in O(M(n)0(n)) operations. 



Proof . Consider the discrete Newton iteration 

\+l = X. -h^(f(x.) -c)/(f(x. + h.) -.f(x.)) . (3.1) 

If hj, = 2""^/^ , ^i"^ " 0[2'^''^) , and the right side of (5.1) is 
evaluated with precision n , then a standard analysis shows that 
x.^ - ^ = 0(2" ) . Since a sufficiently good starting approximation x 
may be found in 0(1) operations, the result follows in the same way 
as in the proof of Lemma 2.2, using the fact that Lemma 2.1 holds with 
M(n) replaced by M(n)0(n) . The assumpfcion C / is only necessary 
because we want to obtain ^ with a relative (not absolute) error 
0(2-") . 

Other methods, e.g. the secant method, may also be used if the 
precision is increased appropriately at each iteration. Q 



I|. Results on Elliptic Integrals . 

In this section we summarize some classical results from elliptic 
integral theory. Most of the results may be found in [1], so proofs 
are omitted. Elliptic integrals of the first and second kinds are 
defined by 

„9 op _T /p 

F(cp,a) = J (1-sin a, sin Q) ' "" d0 (i^.l) 



and 



^ ■ p 9 1 1'? 

E(cp,a) = r (1 - sin a sin 9) / d9 , (1+.2) 



respectively. For our purposes we may assume that a and cp are 
in [0, jt/2] . The complete elliptic integrals, F(rt/2 , cu) and 
E(jt/2,a) , are simply written as F(q;) and E(a) , respectively. 



Legendre's Relation . 

We need the identity 

E(a)F(rt/2-a) +E(jr/2-a)F(ci:) -F(a)F(jr/2-Q;) = jt/2 , (h.^) 

and, in particular, the special case 

2E(jtA)F(rt/^) -(F(V^))^ = ^2 . (^.i^) 



Small Angle Approximation * 

From (i»-.l) it is clear that 

F(q),Q:) = cp+o(a^) {h.^) 

as cc -♦ . 

Large Angle Approximation . 
From (i<-.l), 

F(cp,Q;) = F((p, jt/2) + o(jt/2 -a)^ , (ij-.6) 

■uniformly for < cp <^^ < V^ j ^^ ^ ~* V^ • Also, we note that 

F((p,jr/2) = log tan(Tc/^ + cp/2) . {k.j) 

Ascending Landen Transformation . 

If < a^ < «^+-L < ^2 , < cp^^^ < cPj_ < Jt/2 , 

2 
sin a. = tan (a._^ /2) , (^-8) 

and 

sin(2cp^_^^-cp^) = sin a^ sin cp^ , {k,^) 

then 

^^Vl'^i+l^ " [(1+ sin ap/2]F(cp^,a^) . (1^.10) 



If s. = sin a. and v. = tan((p./2) , then (^.8) gives 

^1.1 = 2sV2/(i,,.) , (,.11) 

and (^.9) gives 

Tj^^^ = Wj/(l+(l+w|)^/2) , (1^.12) 

where 

w^ = ^^ Vl = (v^L+W2)/(l-v^W2) , (J+.I5) 

w^ = tan((p.^^-cp^/2) = w^/ (1+ (1 - w^)^/^) , (If.lH) 

and 

w^ = sin(2cp^^^-cpi) = 2s^v^/(l+v^) . (U.I5) 

Arithmetic -geometric Mean Iteration . 

From the ascending Landen transformation it is possible to derive 
the arithmetic -geometric meem iteration of Gauss: if a = 1 ^ 
h^ = cos a > , 

Vl = (^i^V/^ ' ^^'^^^ 

and 

then 

lim a. = jt/ [2F(a)] . (i^-.l8) 

i -♦<» 



Also, if c = sin a and 



1+1 1 1+1 

then 



(it. 19) 



i=0 



E(a)/F(a) = 1 - L 2^ c7 . (it. 20) 



An Infinite Product / 

Let s. 1 a. and b. be as above, with a = jr/2 -OL^ y so 

1 ' 1 1 ' ' ' 

Sq = b /a . From (1^.11), (i^-.l6) pjid (ii-.J-T), it follows that 
s. = b./a. for all i >0 . Thus, (l+s.)/2 = a.,_/a. , and 

Jf [(l+s^)/2] = lim a^ = rt/ [2F(jt/2-aQ)] (i^.2l) 

i=0 i-*oo 

follows from (il-.l8) . (Another connection between (ij-. 11) and the 

2 2 1/2 
arithmetic -geometric mean iteration is evident if s = (1-b /a ) ' 

Ass-uming (^.11) holds for i < , it follows that 

2 2 1/2 

s . = (l-b./a.) ' for all i > . This may be used to deduce 

(i^-.l8) from (I+.IO).) 



5 • Evaluation of tc ., 

-1/2 
Let a_ = 1 , b^ = c^ = 2 ' . A = lim a. , and T = lim t. , 
'00 -^ . 1 . 1 

i -»oo 1 -♦CO 

where a. , b. and c. are defined by (i4-.l6), (^-l?) a-nd (H.I9) for 

^ n 1 P 
i > 1 , and t. = 1/2 - 1j 2^ c. . From {h.h), (U.I8) and (1+.20), 

^ (0=0 ^ 

we have 

n = A^/T . (5.1) 

Since a. > b^ > for all i >0 ^ and c' ^ = a. -a.^T = a..^ -b. , 
1 — ' 1+1 1 1+1 1+1 1 ' 

(1^.17) gives t.^3_ = [(Vl*°i+l)(Vl-°i+l)l^^^ = Vl-°^=i+l) ' 

2 
so c. p = 0(c ) . Thus, the process converges with order at 

least 2 , and log_ n + 0(l) iteration^ suffice to give an error 

0(2" ) in the estimate of (5»1)- A more detailed analysis shows 



2 2 

that a . . - /t . < jt < a . /"t . for all i > , and also 
1+1' 1 r 1 — ' 

2 i P i + ^P"i 

a./t. - Jt -^ 8j( exp(-2 Jt) and jt - a./t._^ ~ 2 ^ jt exp(-2 jt) 
as i -♦ CO . The speed of convergence is illustrated in Table 1. 



Table 1: Convergence of approxjjnations to n . 



i Tc - a . , T /t . a . /"t . - 3t 

i+1' 1 i' 1 



2.5'-l 8.6'-l 

1 1.0'-3 U.6'-2 

2 7.1^'-9 8.8'-5 
5 1.8'-19 5.1'-10 
ii- 5.5'-lH 3.7'-21 



From the discussion above, it is clear that the following 
algorithm, given in pseudo -Algol, evaluates tt to precision n . 

Algorithm for jt . 

A - 1; B - 2""'-/^; T - l/k; X - 1; 

while A-B > 2"" ^ ' 

1/2 
begin Y - A; A *- (A+B)/2; B - (BY) ^ ; 

T - T - Aa-y)^; X «- 2X 

.^ 

end ; 

2 2 

return A /T [or, better, (A+B) /{hi)]. 

Since logp n + 0(l) iterations are needed, it is necessary to 
work with precision n+0(log log(n)) , even though the algorithm is 
n^umerically stable in the conventional sense. From Lemmas 2.2 to 2.^4-, 
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each iteration requires 0(M(n)) operations, so jt may be evaluated 
to precision n in 0(M(n) log(n)) operations. This is asyinptotical].y 
faster than the usual 0(n ) methods [8,15] if a fast multiplication 
algorithm is used. A high -precis ion computation of jt by a similar 
algorithm is described in [5]' Note that, because the arithmetic- 
geometric mean iteration is not self -correcting, we can not obtain a 
bound 0(M(n)) in the same way as for the evaluation of reciprocals 
and square roots by Newton's method. 



6. Evaluation of exp(x) and log(x) . 

1/2 
Suppose 6 > fixed, and m e [6 , 1-&] . If sin a = m ' , 

we may evaluate F(q;^) to precision n in 0(M(n) log(n)) operations, 

using (^.18) and the arithmetic -geometric mean iteration, as for the 

special case F(it/h) described in Section 5. (When using (k.lQ) we 

need « , which may be evaluated as described above.) Applying the 

ascending Landen transformation (i4-.8-10) with i = 0,1, ...,k-l and 

cp = n/2 gives 

F(cp^,q;^) = I ft [(1+ sina.)/2] iF(aQ) . (6.1) 

1/2 1/2 
Since s^ = sin cc^ = m '^ > 5 ^ > , it follows from (^^-.11) that 

s. -♦ 1 as i -♦ 00 . In fact, if s. = 1-e. , then 
1 11 

^i+1 = ^-'i+l = l-2(l-e.)^/^/ (2-e.) = ^Q + 9(e^) , so s. -> 1 
with order 2 . Thus, after k '-' log^ n iterations we have 
e^ = 0(2"^) , so jt/2-Q:^ = 0(2'^^/^) and, from {k,6) and (h.j), 
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^(V\) = log tan(jr/U + ^^/2) + 0(2""'') . (6.2) 



-rix 



As Sliming k > , the error is uniformly 0(2 ) for all 
m e [& , 1-6] , since ^t. < ^t_ < V^ • 
Define the functions 



U(m) =i Jf [(1 + sina^)/2] |f(o;q) (6.5) 

and 

T(m) = tan(it/l^ + cp^2) , (6.1i) 

where cp = lim cp. . Since s. -* 1 with order 2 , the infinite 

ca .11 
1 -»oo 

product in (6.5) is convergent, and U(m) is analytic for all 

me (0,1) . Taking the limit in (6.1) and (6.2) as n (and hence k ) 

tends to » , we have the fundamental identity 

U(m) = log T(m) . (6.5) 

Using (^-11) to (it-. 15)^ 'we can evaluate 

U(m) =1 ff [(l+s^)/2] >F(q;q) + 0(2"'') and 

T(m) = (1+v )/(l-v ) + 0(2"^^) , to precision n , in 0(M(n) log(n)) 
operations. The algorithms are given helow in pseudo -Algol. 
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Algorithm for U(m) 



.-m) ' ; 
while A-B > 2"^/^ do 



A ^ 1; B - (1-m) 
,-n/2 



begin C ^ (A+B)/2; B ^ (AB)"^^^^ A - C end ; 

A ^ it/(A+B); S - m-^/^; 
rn/2 



while 1-S > 2 / do 



begin A - A(l+S)/2; S - 2S'^/^/(l+S) end ; 
return A(l+S)/2. 



Algorithm for T(ni) . 
V - 1; S *- m-"-/^; 



while 1-S > 2""^ do 



begin W - 2SV/(1+V^); 

' W -W/(1+(1-.W^)^/^); 
W ^ (V+W)/(1-VW); 
V -V/(l+(l+W^)^/^); 
S - 2S"^/^/(l+S) 
end; 
return (l+v)/(l-V) . 

Properties of U^m) and T(m) . 

From (i^■.21) and (6.5)^ 

U(m) = (V2)F(aQ)/F(V2 - a^ , {6.6) 

1/2 
where sin a^ = m ' as before. Both F(a ) and F(7t/2 - a ) may- 
be evaluated by the arithmetic -geometric mean iteration, which leads 
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to a slightly more efficient algorithm for Ij(ni) than the one 
above, because the division by (1+S) in the final "-while" loop is 
avoided. From (6.5) and (6.6), we have the special cases 
U(l/2) = jc/2 and T(l/2) = e^/^ . Also, {6.6) gives 

U(m)u(l-m) = it^/h , (6.7) 

for all m€ (0,1) . 

Althoiigh we shall avoid using values of m near or 1 , it 
is interesting to obtain asymptotic expressions for U(m) and T(m) 
as m -* or 1 . From the algorithm for T(m) , 

as e -» . Thus, from (6.5), 

U(l-e) = L(e) - z/k + O(e^) , 

1/2 
where L(e) = log(i^/e ' ) . Using (6.7), this gives 

U(e) = rt^/[i^L(e)] + 0(e/L^) , 
and hence 

T(e) = exp(ff^/[i^L(e)]) + 0(e/L^) . 

Some values of U(m) and T(m) are given in Table 2. 



1\ 



Table 2: The functions U(m) and T(m) 
m U(m) T(m) 



0.01 


. 6695 


.1.9529 


0.05 


0.8595 


2.5615 


0.10 


0.9821I- 


2.6710 


0.20 


l.i5i^9 


5.1758 


0.50 


1.2972 


5.6591 


0.1^0 


l.ii-522 


I1.I878 


0.50 


1.5708 


1^.8105 


0.60 


1.7228 


5.600^^ 


0.70 - 


1.9021 


6.6999 


0.80 


2.156if 


8.1^688 


0.90 


2.5115 


12.5255 


0.95 


2.87li^ 


17.6617 


0.99 


5.686U 


59.8997 



i i 
Evaluation of exp(x) . 

To evaluate exp(x) to precision n ^ we first use identities 

2 
such as exp(2x) = (exp(x)) and exp(-x) = l/exp(x) to reduce the 

j 

argument to a suitable domain, say 1 < x < 2 (see below) . We then 

/ ' 
solve the nonlinear equation j 

U(m) = X , (6.8) 

obtaining m to precision n , by a method such as the one described 

in Section 5. From Lemma 5.1, with 0(n) = log(n) , this may be done 

in 0(M(n) log(n)) operations. Finally, we evaluate T(m) to 

precision n , again using 0(M(n) log(n)) operations. From (6.5) 

and (6.8), T(m) = exp(x) , so we have computed exp(x) to precision n 

Any preliminary transformations may now be undone. 
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Evaluation of log(x) . 

Since we can evaluate exp(x) to precision n in 0(M(n) log(n)) 
operations, Lemma 5*1 shows that we can also evalute log(x) in 
0(M(n) log(n)) operations, "by solving the equation exp(y) = x to 
the desired accuracy. A more direct method is to solve T(m) = x 
(after suitable domain reduction), and then evaluate U(m) . 

Further Details . 

If X€[l,2] then the solution m of (6.8) lies in (0.10,0-75) , 
and it may be verified that the secant method, applied to (6.8), 
converges if the starting approximations are m =0.2 and m.. =0.7 . 
If desired, the discrete Newton method or some other locally convergent 
method may be used after a few iteration;^ of the secant method have 
given a good approximation to m . 

Similarly, if x e [5,9] , the solution of T(m) = x lies in 
(0.16, 0.85),;^ and the secant method converges if m,. =0.2 and 
m^ = 0.8 . 

If X = 1+e where e is small, and for domain reduction the 
relation 

log(x) = log(\x) -log(x) (6.9) 

is used, for some ?^€ (5^9) j then log(\x) and log(x.) may be 
evaluated as above, but cancellation in (6.9) will cause some loss 
of precision in the computed value of log(x) . If |e| > 2 , it 
is sufficient to evaluate log(Xx) and log(x) to precision 2n , 
for at most n bits are lost through cancellation in (6.9)* On the 
other hand, there is no difficulty if 1 ^ 1 < 2 , for then 
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log(l+e) = e(H-0(2" )) . When evaluating exp(x) , a similar loss 
of precision never occirrs, and it is sufficient to work with precision 
n+ O(log log(n)) , as in the evaluation qf n (see Section 5)' To 
summarize, we have proved: 

Theorem 6.1 . If -» < a < b < oo , then 0(M(n) log(n)) operations 
suffice to evaluate exp(x) to precision n , uniformly for all 
floating-point numbers X€[a,b] , as n -» «> . Similarly for log(x) 
if a > . 



7. Evaluation of Trigonometric Functions . 

-n/2 
Suppose 5 >0 fixed, and X€[6,l] . Let s„ = sin a^ = 2 ' 

2 1/2 
and V = tan((p /2) = x/(l+ (1+x ) ' ) , so tan cp = x . Applying 

the ascending Landen transformation, as for (6.1), gives 

F(q)^.Q^) = j ft [(l+s^)/2]U((P(^,aQ) . (T-l) , 

L i=0 J 

Also, from (l|-.5) and the choice of s , 

'I ^ 

F((PQ,aQ) = artaji(x) + 0(2"^^) . (7.2) 

^ { 

1/2 
From (l|-.ll), S..T > s/ , so there is some j < log^ n + 0(1) 

such that s. €[1/^4-, ^4^/5] . Since s. -♦1 with order 2 , there is 

some k < 2 logp n + 0(1) such that 1-s, =0(2 ) . From (^.6) and 

iS'l), F(^ir^\) = ^^S tan(rt/^ + 9^72) + 0(2"^) . Thus, from (7'1) 

and (7.2), 

artan(x) = / fy [2/(l+s^)] I log taji(V^ + ^^'^) + 0(2"'') . (7.5) 
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If we evaluate tan(jt/^ + ^-^./S) as above, and use the algorithm of 
Section 6 to evaluate the logarithm in {7 -3) y we have artan(x) to 
precision n in 0(M(n) log(n)) operations. The algorithm may he 
written as follows. 

Algorithm for artan(x) , xe[6,l] . 

S - 2"^/^; V - x/(l+ (1+x^)^/^) ; Q - 1; 
while 1-S > 2"^ do 

begin Q - 2Q/(1+S); 

W *- 2SV/(1+V^); 

W -W/(l+(l-W^)^/^); 

W - (V+W)/(l-W); 

V -W/(l+ (1+W^)-^/^); 

S ^ 28-^/^/(1+8) 
end ; 
return Q log((l+V)/(l-V)) . 

k 
After k iterations, Q < 2 , so at most 2 log^ n + 0(l) bits 

of precision are lost because V is small. Thus, it is sufficient to 

work with precision n+0(log(n)) , and Lemma 2.ii- justifies our claim 

that 0(M(n) log(n)) operations are sufficient to obtain artan(x) 

I s 

* \ 

to precision n . 

If X is small, we may use the saiqe idea as that described 
above for evaluating log(l+e) : work with precision 5n/2 + O(log(n)) 
if X >2"^/^ , and use artan(x) =x(l+0(2"^)) if < x < 2""^/^ . 
(Actually, it is not necessary to increase the working precision if 
log((l+V)/(l-V)) is evaluated carefully.) 
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Using the identity artan(x) = jt/2 ~ artan(l/x) (x > O) , we 
can extend the domain to [0,oo) . Also, since artan(-x) = -artan(x) , 
there is no difficulty with negative x . To summarize, we have 
proved the following theorem. 

Theorem "J .1 , 0(M(n) log(n)) operations suffice to evaluate artan(x) 
to precision n , unifonnly for all floating-point numbers x , as 
n -* 00 . 

Suppose 9 € [5 , jr/2 -8] . From Lemma 5«1 and Theorem 1[ *1, we 
can solve the equation artan(x) = 9/2 to precision n in 
0(M(n) log(n)) operations, and thus evaluate x = tan(9/2) . Now 
sin 9 = 2x/(l+x ) , cos 9 = (1-x )/(l+x ) , etc., may easily be 
evaluated. For arguments outside [8, jt/2-8] , domain reduction 
techniques like those above may be used. Difficulties occur near 
certain integer multiples of jt/2 , but these may be overcome (at 
least for the usual floating-point number representations) by increasing 
the working precision. We state the following theorem for sin(x) , 
but similar results hold for the other trigonometric functions (and 
also, of course, for the elliptic integrals and their inverse functions). 

Theorem 7.2 . If [a,b] c (-Jt, jt) , then 0(M(n) log(n)) operations 
suffice to evaluate sin(x) to precision n , uniformly for all 
floating-point numbers xe[a,b],asn-»co. 
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Q* Asymptotic Constants . 

So far we have been concerned -with order of magnitude results. 
In this section we give upper bounds on the constants K such that 
w(n) < (K+o(l))M(n) log n , where w(n) is the numher of operations 
required to evaluate jt , exp(x) j etc., to precision n . The 
following two assumptions will be made. 

1. For all 7 > and e > , the inequality M(7n) < (7+e)M(n) 
holds for sufficiently large n . 

2. The nimiber of operations required for floating-point addition, 
conversion between representations of different precision (at 
most n ), and multiplication or division of floating-point 
numbers by small integers, is o(M(n)) as n -^ oo . 

These assumptions certainly hold if a standard floating-point 
representation is used, and the multiplication algorithm has 
M(n) ^n(log(n)) (log log(n))'^ for some a > , provided p >0 
if a = . 

The following result is proved in [5]« The algorithms used are 
similar to those of Section 2, but slightly more efficient. 

Theorem 8.1 . Precision n division of floating-point numbers may 

be performed in (^1-+ o(l))M(n) operations as n -» «> , and square roots 

may be evaluated in (11/ 2 + o(l))M(n) operations. 

Using Theorem 8.1 and algorithms related to those of Sections 5 
to 7^ "tlie following result is proved in [1^4-]. 
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Theorem 8»2 . jr may be evaluated to precision n in 
(15/2 + o(l))M(n) log n operations as n -* ^ . If jr and log 2 
are precomputed, the elementary function f (x) can be evaluated to 
precision n in (K+ o(l))M(n) log n operations, where 



K = 



15 if f (x) = log(x) or exp(x) , 

5^ if f(x) = artan(x) , sin(x) , etc.. 



and X is a floating-point number in an interval on which f(x) is 
defined and bounded away from and <» . 

For purposes of comparison, note that evaluation of log(l+x) 
or log((l+x)/(l-x)) by the usual series expansion requires 
(c+ o(l))M(n)n operations, where c is a constant of order unity 
(depending on the range of x and the precise method used) . Since 
15 logp n < n for n > 85 ^ the 0(M(n) log(n)) method for log(x) 
could be faster than the 0(M(n)n) method for n greater than a 
few hundred. 
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